#include "cuda_runtime.h"
#include "cublas_v2.h"
#include <stdio.h>
#define MACH_EPS 2.220446049250313e-016
#define IDF2(i, j, K) ((j)*(K) + (i))
#define min(X,Y) (((X)<(Y))?(X):(Y))
#define max(X,Y) (((X)>(Y))?(X):(Y))
#define abs(X)   (((X)>=0)?(X):-(X))
#define sign(X,Y) (((Y)<0)?-abs(X):abs(X))

#define checkCudaErrors(err)           __checkCudaErrors (err, __FILE__, __LINE__)
#define getLastCudaError(msg)      __getLastCudaError (msg, __FILE__, __LINE__)


typedef enum
{
	DEPT150_SUCCESS = 0,
	
	DEPT150_MACHINE_SINGULAR = -1,

	DEPT150_POSSIBLE_SINGULAR = -2,

	DEPT150_NOT_POSITIVE_DEFINED = -4,

	DEPT150_CUBLAS_ERROR = 1000,

	DEPT150_MKL_ERROR = 2000
} dept150_error_t;

typedef struct gpu150_device_matrix {
	int order;//order of the matrix
	int block_size;//block size
	double** gpu_matrix_blocks;//array of pointers to GPU memory space where matrix is stored
	int gpu_n; // number of gpu used for storing matrix
} gpu_based_matrix_t;

inline void __checkCudaErrors( cudaError err, const char *file, const int line )
{
	if( cudaSuccess != err) {
		fprintf(stderr, "%s(%i) : CUDA Runtime API error %d: %s.\n",
			file, line, (int)err, cudaGetErrorString( err ) );
		exit(-1);
	}
}

inline void __getLastCudaError( const char *errorMessage, const char *file, const int line )
{
	cudaError_t err = cudaGetLastError();
	if( cudaSuccess != err) {
		fprintf(stderr, "%s(%i) : getLastCudaError() CUDA error : %s : (%d) %s.\n",
			file, line, errorMessage, (int)err, cudaGetErrorString( err ) );
		exit(-1);
	}
}

int gpuDeviceInit(int devID);
void print_matrix(int rows, int cols, double * matrix, int ld);
void printGPUMatrix(gpu_based_matrix_t gpu_matrix);

int copyMatrixToDevices(double *matrix, gpu_based_matrix_t *gpu_matrix, cublasHandle_t* cublasHandles );
int copyMatrixFromDevices(double *matrix, gpu_based_matrix_t gpu_matrix );

double * getA15(int order);//cond 4*n*n
double * getA16(int order);//cond 2*n*(n+1)

//calculate Matrix norm by columns
cublasStatus_t getMatrixNorm( gpu_based_matrix_t gpu_matrix, cublasHandle_t* cublasHandles, double *norm );
//get matrix cond
cublasStatus_t findMatrixCond(double *host_matrix, gpu_based_matrix_t gpu_matrix, cublasHandle_t* cublasHandles, double matrixNorm, double * cond, bool isSymmetric, bool isUnitDiagonal);

int getGPUCoordinate(gpu_based_matrix_t gpu_matrix, int row, int col);

int initCublasHandles(cublasHandle_t* handle, int gpu_n);
int destroyCublasHandles(cublasHandle_t* handle, int gpu_n);

//int initCudaStreams(cudaStream_t* stream, int gpu_n);

void solveLTempEqualsB(double *host_matrix, double *solution,  int order, int *ipiv, bool isUnitDiagonal);
void solveLTempEqualsB(double *host_matrix, double *solution,  int order, bool isUnitDiagonal);
void solveUXEqualsTemp(double *host_matrix, double *solution,  int order, bool isSymmetric);

cublasStatus_t solveLowerTriangularSystem(gpu_based_matrix_t gpu_matrix, cublasHandle_t *handle, cudaStream_t* stream, double *solution, int *ipiv, bool isUnitDiagonal);
cublasStatus_t solveUpperTriangularSystem(gpu_based_matrix_t gpu_matrix, cublasHandle_t *handle, cudaStream_t* stream, double *solution, bool isSymmetric);

